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SUMMARY | 


The use of computers to optimize free parameters of a system has become 
relatively widespread in many areas of engineering. Parameter optimization 
codes have been written for that purpose, and make it possible for a design 
engineer, once he has developed the mathematics of a system, to optimize its 
parameters according to some criteria. But of equal interest to the design 
piyaneer is the sensitivity of the optimized criteria to departure of the 
parameters away from their optimal value. The purpose of this thesis is to 
show ways in which a parameter optimization code may be augmented to yield 
Such sensitivity information. 

A Fletcher and Powell version of Davidon's variable metric optimization 
search technique was employed to optimize multi-parameter functions. Their 
method is useful in that it computes the inverse Hessian matrix (or matrix of 
second derivatives), which completely describes the curvature of the function 
at the optimum. Equations were developed so that the sensitivity could be 
expressed in a meaningful output format. This was made possible through the 
use of matrix inversion and eigenvalue analysis subroutines which were ob- 
tained from the scientific subroutine library of the IBM 360 91 and used in 
conjunction with a digital computer code employing the Fletcher and Powell 
technique. Equality constrained optimization problems were also considered 
by employing the penalty factor method proposed by Courant and used by Kelley. 
Equations analogous to the use of Lagrangian Multipliers were used to deter- 
mine the cost of the equality constraint. 

Example problems are offered showing the optimal solutions, sensitivity 
data at the optimum, and interpretation of that data. The well known 
Rosenbrock function was used to exhibit the accuracy of the methods employed. 
A typical engineering problem was solved involving the sensitivity of an 
optimal nuclear rocket engine used to inject a payload onto an interplan- 
etary trajectory. The results indicated that the thermal power of the re- 
actor and the ratio of length to diameter of the core could be varied con- 
siderably from their optimal values with little cost. The power density 
however was relatively fixed for optimal operation. 
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Te, INTRODUCTION 


The use of computers to optimize free parameters of a system has be- 
come relatively widespread in many areas of engineering. Parameter opti- 
mization codes have been written for that purpose, and make it possible 
for a design engineer, once he has developed the mathematics of a system, 
to optimize its parameters according to some criteria. But of equal in- 
terest to the design engineer is the sensitivity of the optimized criteria 
to departure of the parameters away from their optimal value. The purpose 
of this thesis is to show ways in which a parameter optimization code may 
be augmented to yield such sensitivity information. The interest in this 
is evident. 

Since the engineer is working with temperatures, pressures, masses, 
etc., the computed optimal solution cannot be implemented exactly. Physical 
parameters are subject to uncontrollable variations and the resulting de- 
parture from the optimum may be quite significant. For space flight 
applications, maximum payloads might be of primary interest for one mission 
while a minimum fuel consumption the criterion for the next. Design require- 
ments for a Venus flyby will certainly differ from those of a Mars lander 
or Jupiter probe, and yet many of the systems must be flexible enough to 
be useful for all three. Thus, the optimized solution must also contain 
significant information about departures from the optimum. This thesis 
addresses that problem; i.e. sensitivity analysis at the optimum. 

Parameter optimization algorithms are many in number and varied in 
application. They differ primarily in their rate of convergence and the 
restrictions imposed on the function under consideration. Nearly all 


require a large number of iterations to achieve a given accuracy in locating 





the optimum, and some procedures may not converge from a poor starting point. 


Because, near the optimum, the second order terms in the Taylor series 
expansion dominate, the only method which will converge quickly for a gen- 
eral function are those which will guarantee to find the optimum of a 
general quadratic function speedily. Fletcher and Powell [1]* have pro- 
duced such a method which was based upon a procedure described by Davidon 
[2]. The method is superior to other techniques which possess quadratic 
convergence in that it makes use of information determined by previous 
iterations and also in that each iteration is quick and simple to carry 
out. The primary justification for its use however, lies in the fact that 
the method yields the necessary information to determine the curvature of 
the objective function at the optimun. 

The method of Fletcher and Powell estimates and continually updates 
the inverse of the Hessian matrix, H, (or matrix of second derivatives) 
during the optimization search so that a close approximation to the true 
value at the optimum is peaches Through the use of the eigenvalues and 


eigenvectors of the Hessian matrix, analysis of the characteristics of the 


objective function in the neighborhood of the optimum is possible. 


As previously mentioned, it is of interest to know not only where the 


optimum is located, but by how much each parameter Xx, may be changed from 


* Figures in parenthesis indicate references listed following the text. 








its optimal value before a significant change, A, in the objective func- 


mlon occurs. The variation of jthe x, may be independent of the other ,., 
al 

or several parameters may be changed simultaneously in a co-ordinated fashion 

away from the optimum. The latter variation might allow for significantly 


large departures from the optimum for a specified A when the function 


pa + von ° = imi * 3 {ft 
en W-eimensione ridge . 
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In this respect, it is intended that this thesis may be used to evalu- 
ate and analyze the optimum of any general function of N independent 
variables in such a way that a complete sensitivity at the optimum is 
clearly presented in a useful output format. It is shown to the user 
which of the specified variables may be changed and the magnitudes of those 


changes before specified degration in the objective function will be 


exceeded. 





II. PROBLEM DEFINITION 


An analysis of the behavior at the optimum of an N-variable function 
is possible where the second derivatives are known. Suppose that Y is 
a real valued function of N variables with continuous first and second 
partials and possesses a relative minimum at Xo then the first deriva- 


tive will vanish and by Taylor series expansion: 


Y-Y = AY == Ax be AX, + Higher Order Terms (1) 


=> 4 
where i is the Hessian matrix (H) of Y at the optimun. 
Ox 


From the Taylor series approximation (1) we find that the gradient 


Ox 


ic a 
AY(X) = ree oo) VY (X) (2) 


and solving for x. yields: 


=} 
X =X- ns 0) VY (X) (3) 
Oo ane 


-1 

2 

Se that if rea) were known, the step to the optimum would be 
Ox 


given by (3). Some optimization algorithms collect information which 





generate an approximation to the inverse Hessian matrix during the search 


—1 


for the optimum. These algorithms provide H So) that’ the probuienvos 


Sotaining HW, in the use of such an algorithm is reduced to the relatively 
fa 2s ok 

simple additional task of inverting the Nx N matrix ot 
ny 2 
ox 


Historically the method is similar to Newton's method which minimizes a 
function Y(x), x on N-vector, by generating a sequence of points (compare 


with equation (3)) 


y OL) 2 yk) _ Ck) 0k) (k) 
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Where is an appropriately chosen scaler constant. Fletcher and 


Powell's version is in fact a quasi-Newton method which uses an initial 
; k)4- 
estimate and computational history to generate an estimate to cul dy : 


rather than performing the computational work of evaluating and inverting 


the matrix at each step. 





III. DAVIDON'S METHOD 


Davidon introduced a variable metric method which was the first 
optimization technique to use past information estimating the inverse 
Hessian matrix. Fletcher and Powell have improved upon the method by 
simplifying the iteration scheme and modifying the criterion for con- 
vergence. Their method, which numerically determines the minimum of 
a function of several variables, combines the best features of the con- 
ventional gradient method and Newton's method, namely the sureness of 
convergence of the former and the quadratic terminal convergence of 
the latter. An excellent exposition of the method, including conver- 
gence proofs, has been given by Fletcher and Powell. For the purpose 
of this thesis however, it will suffice here to state the algorithm and 
point out the usefulness of its main features as was done by Kelley [3]. 

Let us denote the free parameters as x. and the function to be 
Minimized as Y(x). It is assumed that Y has continuous first and 
second partial derivatives with respect to xX. Any starting point x” 
is chosen according to some criteria*. At the starting point ae the 
gradient vector well as Y itself, is evaluated. A change is 


, 0 , 
then made in xX according to: 


AX = -a HY (4) 


x 


* The freedom in the choice of the starting point depends on Y. If Y 
is well behaved, this choice is free. In other cases, the starting 
point must be close to the optimum to assure convergence. 





where H™~ isa positive definite, symmetric matrix, defining a metric 
gm the X-hyperspace. Ets initial selection iswotherwise arbi trary ie 
general purposes the unit matrix may be used, but if the parameters differ 
by orders of magnitude it is convenient to either scale them or to estimate 
the accuracy with which they are to be determined. a& > 0 is a step size 
parameter. 

In Davidon's method, the one-dimensional minimum of F vs. a is 
obtained along the vector originating in x° in the direction H™ 


At the new xX, the gradient vector uo is again evaluated. The a 


matrix is updated according to 


Ho? + AH? =H? +A+B (5) 
where 
oly 
4 = AX_AK 
AX AY 


and generates the inverse of H in N _ steps for the quadratic function 


and where 


is intended to cancel the initial assumption for HO {4 ]. The procedure 


. , ~1 
is begun again with the new values of X, Yo» Hoe 





It is shown in [1] that H remains positive definite and that, as 
X approaches the minimizing point, H™? approaches oe evaluated at 
the minimum. For quadratic Y, in N-dimensional space, the minimum is 
obtained in, at most, N steps (within round-off error): the method is 
quadratically convergent. 

For more general functions having the smoothness properties assumed, 
a Taylor expansion through quadratic terms provides a good representation 
of the function in some neighborhood of the minimum. With H™ converged, 


mc minimum Of Y vs. @ then will occur for a = 1, 


KOS =o Gf eey 
x 


will approach the value given by Newton's method, namely = ey Ye 





IV. SENSITIVITY OF THE PERFORMANCE FUNCTION AT THE OPTIMUM 


i 


Assuming that the optimum and H ~~ at the optimum have been deter- 


mined and that H? has been inverted, then all the necessary information 
has been collected for a detailed analysis of the sensitivity of the per- 
formance function to changes in the performance parameters. 


bec the criterion function be 


ema assume for simplicity that the origin, i.e. 


is an analytic optimum of Y. Around the optimum 


Y -5 x "ax (6) 
where 
! d2y 
H = - la | 
Ox. OX 1 


t= the Hessian matrix of Y at 0. 





Independent Variation (one at a time) 


If one parameter at a time is allowed to depart from the optimua, 


the function Y is 


where a. 1s the corresponding diagonal term of H. For an assigned 


change A in Y we find 


+ Ax, = : | 


ear \ 
2AY 


re Ae e 
aa 8 








(7) 


where +i Ax. is the allowable change in Xx. for the previously assigned 


acceptable variation in Y. 


Simultaneous Variation 

If one allows several parameters to be changed in a co-ordinated fashion 
away from the optimum, departures from O for a specified AY may be 
significantly larger than those shown in the one-parameter-at-a-time case 
To illustrate this, consider the case in a 2-dimensional parametric space 
where the function y presents the characteristics of a ridge as illustrated 
imefigure 1. 

The Y =lA line intersects the x, and x axes at distances which 


1 Z 


mee tar less than the values of x4 and x, at the end of thew y i= 0) 


contour. (Five times less in the case of Figure l). 
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A Skew Ridge Illustrating Advantage of Simultaneous Variation 
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Eigenvalues and Eigenvectors of al 


To be able to analyze the characteristics of the Y surface in the 
neighborhood of 0, it is convenient to make use of the eigenvalues and 
eigenvectors of H at that point. 

It is assumed here that the final use of this analysis will be in 
computer programs, and that matrices eigenvalue and eigenvector search 
smoroutines are G@ither available or easily implementable. To that respecr. 
we note that H is a symmetric, positive definite matrix, and that know- 
ledge of this fact simplifies the implementation of such subroutines. 

Let d. and Vv, be respectively the eigenvalues and eigenvectors 


mau f 2], i.e. 


det (H - AI) = 0 


Since H is real and symmetric the eigenvalues are real and the 


eigenvectors are orthogonal and may be expressed in orthonormal form: 


Vv ~V~. = 1 (normal) 
i i 
ec aes = 
= : Ve =Q i # 4 (orthogonal) 


If the N-parameters are allowed to be changed in the direction of the 


Ji eigenvector, i1.¢e. AX =k. V., then the degradation in the perfor- 
5 et 


mance function AY is: 


heen 





AY Li we V 
ae c. i | Kk. Le (8) 
pit 
= ie a oe ae 
HV. =A, V. and V.*V. = 1 
i i a i 
so that 
ae 2 
Me 5 K, d. (9) 


Emaeror an allowable variation A in Y 


_ PX 
k. = e CLO) 


: Ben 
k. represents how far in the direction of the normalized i eigenvector 
one may travel on the response surface before degrading the performance 


function by AY. It is observed that k is a maximum for A win? 


SO that the direction of least 


sensitivity to changes in the 

parameters, Be: is given by 

the eigenvector with A=A.. 
min 


In other words the length of the AX vector = is maximized for 





a given’ AY in the performance function. 
Consider a hypothetical 2-variable optimization problem where the eigen- 


values of H at 0 assume the values of 1 and 10, i.e. 


=o= 





and the associated orthonormal eigenvectors are: 
v, =e 0) and V5 = (On, 


Bespectively. The contour of the objective function at the optimum wall 
then assume the shape of Figure 2. 
Observe that for a given Ay, the distance from the optimum in the 
direction of We is considerably greater than in the direction of V 
[ 52 


—£ 
‘a 


> 
imaetact, it follows directly from (10) that the distance is 


\ 
‘ 


ereater. 


It is evident that the relative length of any existing ridges in the 
response surface will be determined by the square root of the ratio of 
any two eigenvalues and that the magnitude of the sensitivity for a 
given A will be determined by the square root of the eigenvalues. 
Let two distinct solutions of optimization problems assume the values: 
Solution 1 (as before) SolMErone2 


Ava = 1 dio = 10 hoy = 0.1 Ano = 1 


= (0,1) ¥, = (1,0) V, = (0,2) 


The response surfaces will then assume the shapes of Figure 3a and Figure 
3b respectively. 

Notice that the ratio of the eigenvalues has remained constant and 
that the shape of the response surface is unchanged. The magnitude of 
the eigenvalues was decreased by a factor of 10 and if drawn to scale the 


response surface would be expanded in all directions by a factor of Y10. 


p= 
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Effect of Eigenvalue Ratio on Response SULLaece 


Figure 2 


~|5- 








Figure 3a 
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Effect of the Magnitude of Eigenvalues on the Shape of the 


Response Surface 


Figure 3 
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Skew Ridges 

The previous examples have pointed out that when the matrix H has 
one eigenvalue much smaller than all other ones a ridge will exist in the 
response surface. A skew ridge exists when one ecigenvalue is much smaller 
than all other ones, and the corresponding eigenvector is not parallel to 


one of the axis in X-space. 


iiyaecrdeé Vsmparavlel tovone 


axis in X-space, it can always 


be removed by fal change in scale 
along that axis. Jt is there- 
fore Lepreseneaurve Or Sine way 


scales are chosen rather than 


a characteristic of the func- 


ELON EOsbe Opel Zed - 


The situation is completely different if the ridge is at an angle to 
the axis because no change in scaling can remove the ridge. Such a ridge 
reflects an interaction between parameters in the way in which they affect 
Ene Criterion function. 

Only when the eigenvalues and eigenvectors of the Hessian matrix are 


known, can such a ridge be discovered; and only then can the characteristics 


of the optimum be determined. The complete search of eigenvalues and eigen- 
vectors of H and the derivation of the resulting sensitivity information 


iS contained in the computer program described in the following section. 


S172 
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V. COMPUTER PROCRAM 


The program described in this thesis makes use of a standard, 
Powell-Fletcher type, paramecter-optimization computer software 
package (6). It contains a multi-dimensional optimization algorithm 


similar to Davidon's variable metric method as was previously described. 


The matrix inversion and eigenvalue analysis subroutines were obtained 
meom che scientific subroutine library of the IBM 360 model 91. A listing 
of these subroutines is presented in the Appendix. A simplified flow 
diagram illustrating the coordinated use of these programs is shown on 


the following page. 


Necessary Input for Sensitivity Output 
The input is the expression relating the objective function Y to 
the N independent variables Xs The expression may take the form of 


a Single equation, e.g. 


Y= £(x,) 


or may comprise any number of subroutines as long as the objective function 
Y and the N independent variables are defined. The allowable departure A 
from the optimum is also a required input and must be specified by the user. 


It may be expressed as a percentage change in Y, e.g. 





Y = f(s) 


oe ~ 
wf 7 
—/ 
eae td 
- -_ a A a — en 


ice = 
POWELL - FLETCHER 
OPTIMIZATION ALGORITHM 
— penis 2". 


| 





x Y i 
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EIGENVALUES, 
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Si 
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| EQUATIONS . 
ic 


AX. (independent) 


gG— 


(AX). (along eigenvector) 


Flow Chart of Computer Programs 


= 10= 





AY = iby ee 
or Y = 99% Y 
opt 
or as a fixed quantity, e.g. 
AY = 10 
Output Format 
Given the required input: 
Y = E(x, ) 


Y=A (a Conetane) 


the sensitivity information is presented in the following table for con- 


venience in analysis. 


TABLE I 
AxvEOr Speci led sy 
: AY x. Independent Simultaneous Variation 
OP opt Variation KM, SS = a ee = 
mn max 
=. oe ee 
x, Ax, Ax, Ax, 
su 
Xo t+Ax. Ax, Ax, 
Xv Ax, Ax. | Ax. 


~2?0- 





Eo 


VI, NUMERICAL EXAMPLE: ROSENBROCK FUNCTION 
In order to illustrate the usefulness and accuracy of sensitivity 
analysis a standard test function is offered as an example. The function 


- Me * = 2 
Y (x 100 (x., xy )* + (1 x,) 


me 
I? 2? 
was proposed by Rosenbrock and is interesting in that it possesses a 


Steep-sided ridge following the curve x * =x, as shown in Figure 4. 


I! 2 
The exact solution of the problem is offered along with the computed 
values so that a comparison can be made and the effectiveness of the 
techniques employed may be evaluated. 


Problem: minimize Y = 100 (x, ~ ca) + (1.0 - oi 


Bact Solution: 


Y = 0.0 
opt 
802 -400 0.5 1.0 | 
ae 10 H = H? = 
-400 200 1 2.005 
my So ee 


where all numbers given are exact. 


Fletcher and Powell's version of Davidon's variable metric technique 


found an optimum solution: 


= 107° ~ 0.0 
opt 
es 1.00007 = 1.0 
x =O eG 
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Rosenbrock's Curved Valley Function 


ye 





Figure 4 





and the approximated inverse Hessian matrix 


025055 1.0051 
ee 
OOS 2.0035 


which represents an error of less than 12. 
The matrix inversion subroutine with the approximated H™ as its 


input resulted in the Hessian matrix: 


S202 =415.5 


= 555 209.0 


ru error of about 4%. This error is more than acceptable for the purposes 


of sensitivity analysis. 


> 


TABLE IIL 


Sensitivity Data of Rosenbrock function: 





hoor Gor opecimacd oy 


F AY x. , Independent Simultaneous Variation 
a opt Variation eel 399 Kae MOST 
min 2 
on” 0.10 1.00007 pe O15 Her ceil) | @ lbe  a2. OL 
/ 
HOOOL7 0309 Ax. = +.633 Ax. = -.006 
| 
Analysis: 


If the response surface were to be constructed from the sensitivity 
data shown with no knowledge of the function under consideration it would 


resemble Figure 5. 
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Predicted Shape of Rosenbrock Function at the Optimum 


Figure 5 
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A comparison with Figure 4 shows that the sensitivity analysis at 
the optimum gives the desired results. The above figure shows that the 
mincLiOn 1s extemely sensitive fo changes along Eheux) axis (Ax, =art iO) Le) 
for Y = 0.1) as well as to changes along the * axis (Ax., = +, 03 eee in 
contrast, it is shown that a skew ridge exists along the direction defined 
by Ax, = .318, Ax, = .633 and that if varied simultaneously x, may 


2 
be changed by 304 and Xo by 60% before Y changes by 0.1. 
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VII. EQUALITY CONSTRAINTS 


In the optimization of multi-variable problems it is often the case 
that only certain combinations of the variables are either meaningful 
or acceptable. The imposed restriction usually assumes the form of an 
equality (or inequality) constraint: 


G(x.) = 0 
ot 


The analytical solution of such a problem can be found by the method 
of Lagrangian Multipliers [7], which seeks values of the parameters for 
which the modified objective function 

F* = F + AG GEES 
is stationary, 1.e. 


F* =F + AG = 0 (125 
x x 


Solving for 2A yields 


rj 


x 
A= - z (13) 
x 


X is meaningful in that it represents the cost of the constraint G. If 


G were relaxed by 1 unit then F could vary by i. 





In the general case however, numerical search methods are unlikely 
to locate all types of stationary points for the modified function using 
Lagrangian Multipliers. A more feasible approach is that of penalty 


Bactors. [8] 


Benalty Factors 


In the use of penalty factors a modified function which incorporates 


ne constraint is defined as 
Fx = F + PG* (14) 


where P is a large positive-valued constant (for minimization). If P 
is chosen large enough then G is forced close’to zero in the search for 
the optimum. At the optimal solution G is equal to some small quantity 
€. If e€ is not within an acceptable distance from G then P is in- 
creased until an optimum i found which satisfies G within e¢. The 
cost of the constraint can be found in a manner identical to the method 


of Lagrangian Multipliers. 


At the optimum of the modified function of (14) 


F* =F + 2PEG =0 (i>) 
x x Xx 
so that 
F 
x 
2Pe=-G =A (16) 
x 


is the cost of the constaint G. 
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Although the idea of a penalty function seems to have been conceived 
some years ago (Courant [1943]), there has been very little computational 
experience with regard to its applications. By examining equation (16) 
it is observed that as € approaches zero, P becomes infinitely large. 
Large values of P however effectively produce a sharp ridge in the con- 
tours of F*, and most search techniques are troubled by the existence 
Seesuch@a ridge. The choice of P is therefore a compromise between 
large values for small violations in G, and smaller values to eliminate 
troublesome ridges in the modified objective function. 

In using Davidon's optimization technique however, it was discovered 
that even if a ridge presents no problem in finding the optimal solution, 
P may not be chosen arbitrarily large. For this situation, € becomes 
so small that changes in the parameters of order e€ produce corresponding 
changes in F which are less than the criterion for convergence in the 
search for the optimum. it € is to be meaningful it must be large enough 
to possess a unique solution. In other words changes in € must be large 
enough to affect the terminal convergence of the optimization search. It 
is therefore necessary to possess some insight into the problem before 
the penalty factor method can be employed. 

We may note that for any given equality-constrained optimization 
problem,. A4 will possess a unique value. Analytically A is found to 
be 2PE€ as € approaches zero and P approaches infinity. Let us 
denote oe as the maximum allowable violation in the equality constraint 
and ein 3 the minimum e€ which will produce a unique optimum within 


the limits of the convergence criterion. € must now satisfy 


ewes € Se 
m™ 


mith ax 


for a meaningful solution. 
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Values for A and P are estimated within one or two orders Gr 


magnitude such that 


BE seas) ae Nae 


If the resulting optimum possesses an allowable e€ then the solution is 


found with the cost of the constraint 
XN = 2Pe (17) 
and the constraint violation: 
gE =G = G (18) 


If €O€E Orweens € . then P must be increased or decreased re- 
max ee gh 


Spectively until an allowable e€ is found. 


Example Problem 


To illustrate the use of the penalty factor method for equality con- 
strained optimization problems the following two examples are offered for 
comparison and analysis. 

m _ 2 2 _ 2 
er Fx, Xo) = F(X» Xo) (x, 3)* + (x, 2) 
and CG, Gr» xo) =x, + xy 4 = 0 
G (x, > x 


g? = ®1 %o 


Zyros 





Constrained minimum 





Figure 6a 





Figure 6b 
Bie 
Meer pEectation Of eCOSE O© EherConstraintm — 


G 
ae 


Figure 6 
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The response surfaces are shown in Figure 6, where the equality 
constrained function has a minimum value at the same point for both 
Beoplems, 1.¢e. xX) =X, = 2.0. Sinee the unrestricted functionercmene 
same and the solution lies at the same point, then the gradient F is 
the same for both problems. The constraint and its gradient G_ are 
different however and consequently the cost of the constraint A will 
have two distinct values. 

Analytically it can be shown that SELUGs =. = 2 for the linear 
constraint x, + x, = 4 and that “F/G, =} = 1 for the hyperbolic 


Senstraint xy X, = 4, The cost of the constraint has been reduced in 
the second case because F is less sensitive to changes in G, at the 


optimum. 


Computer Results: 


For F* = F+P rn and F*¥ = F +P G,? 
and P =100 P = 100 
F* = 1.980 F* = 1.994 
x, = 2.005 x, = 2.0009 
x, = 2.005 X, = 2.0020 
ence = C= 200 —E = .0057 
A = 1.988 X = 1.141 


In the second case, a minimum was found which was very close to the 
constraint, e.g. € = .0057. As aresult A =1.141, a difference of 
.141 from the analytical value. This error may be explained by examining 


the magnitude of ¢€. An e€ of .0057 is very small for the problem under 
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consideration. Changes in the variables of order € will not affect the 


optimal solution, so that € is actually less than as previously 


“MIN 
defined. 
A second optimization was performed with P reduced to 50 in order 


to increase €. (Remember that it has been shown that 2Pe is constant). 


The results were: 


les SEAS 8) 
x, = BSUS s 
X= 999 
2 = lea: 
Xr = 1.07 


It is observed that decreasing P resulted in a more meaningful value 


for € and consequently a closer approximation to the const constraint, A. 


=2 = 
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VIII. SENSITIVITY ANALYSIS OF A NUCLEAR ROCKET ENGINE 


As was previously mentioned, the primary purpose of a sensitivity 
analysis is to aid the design engineer in constructing a system which 
ful operate in an Optimal fashion, despite variations in the coneroiiane 
parameters. To that effect, an illustrative example is presented whereby 
the engine design parameters of a nuclear rocket are optimized in order 
to achieve a maximum payload at a specified hyperbolic velocity. 

A set of mathematical models of the elements of nuclear rocket 
engines, suitable for detailed systems analysis, has been developed [9] 
which constitutes the basis for a digital computer program called 
NUROC/SAC (Nuclear Rocket Systems Analysis Code). The Code has been used 
to describe a number of existing engines and the results obtained weve 
found to be accurate [10]. 

ESCAPE [11] is another computer code which calculates geocentric 
(or planetocentric) tangential thrust escape trajectories and which may 


be used in conjunction with NUROC/SAC. 


Input Parameters 


From a design viewpoint the most important input parameters to 


NUROC/SAC are: 


Q thermal power of the reactor, watts 
D diameter of the reactor core, meters 
L length cf core, meters 

L/D ratio of core length to diameter 


maximum allowable core material temperature, °K 
cmax 





pe nozzle chamber stagnation pressure, n/cm* 


NAR NOZzbe wane Catlo.: a5 /A 


exit fmrooul 


The output format of NUROC/SAC consists of a summary of operating variables 
including the input design parameters, defining entirely the characteristics 
of a specific nuclear rocket engine. The most important of these and the 


ones which will be used as inputs into ESCAPE are: 


e 


m total engine propellant mass flow 
F total engine thrust 
mi, total engine mass 


Additional inputs to ESCAPE which must be specified are: 


mo initial mass in earth orbit 

Hy as orbital altitude 

Vue final hyperbolic excess velocity specified 
Ke tankage to propellant mass ratio 


The output of ESCAPE may assume a variety of formats but the payload 
delivered at the specified hyperbolic velocity is the payoff in the opti- 
mization problems which follow. The payload is defined as the initial mass 
minus the engine mass, the tankage mass, and the propellant mass required 


to reach Vie’ 


A previous study [12] has indicated that within the range described by 
technological constraints the payload performance of the nuclear rocket 


will increase monotonically with the maximum allowable core material 


Sie 
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a rrrvenmnyss «= ANALYSIS 








ENGINE PARAMETERS 
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temperature and the nozzle area ratio, and will decrease monotonically 
with the chamber stagnation pressure. Thus choice of these parameters 
is limited to technologically realizable values. Typical values presently 


ii use are: 


= 2000 7k 
cmax 
P. = 300 n/cm? 
NAR = LOOR 


The power, power density, diameter, and length which parameterize the 
core geometry may be varied in a constrained but optimal fashion to describe 
an engine which will inject the maximum payload into a specified inter- 
planetary trajectory for a given initial mass in earth orbit. The problems 
which follow will consider an initial mass of 100,000 kilograms in an earth 
orbit of 500 kilometers and will optimize the core geometry to inject the 
maximum payload into a trafectory described by a hyperbolic excess velocity 


of 10,000 m/sec. 


Optimization of (Core Geometry 


ime is interesting to mote that 1f any three of Ene four eneine 
parameters, power (Q), power density (WV), diameter (D) and length (L), 
are specified then the fourth is automatically determined. Since the 
length and diameter describe the volume, the power density can be ex- 


pressed as a function of Q, L, D, i-e., 
Z. 
We e407 Wel 


Thus the optimization of core geometry is reduced to a 3-dimensional 


problen. 
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Maximize Payload: Free Parameters Q, w, L/D 
The following example maximizes the payload delivered by a nuclear 
rocket with free parameters: power (Q), power density (W), and ratio 


of length to diameter (L/D). 


Beortineg Point; 


Q = x, = 2000. MWt 
yp = X= 6.0 x 10? W/m? 
L/D = = 4.0 
/ X. 
TABLE (iit 


Sensitivity Data at the Optimum 





Ax for Specified Ay 


1a : Sy Ga apis X - Independent Simultaneous Variation 
P oP Variation A => en = Se X = 6150 
31,630 Kg.| 316. Kg | Q=1832 MWt | AQ=+333. Gor -28. 23% 
W=8.84x10" | Av=4+3.034 | 3.88 -.062 g -.017 , 
W/m? (x107) (x107} (x10°) (x10 
L/D=4.80 AL/D=+.57 ood 677 «Is )S) 


Analysis: 

When both wy and L/D are free parameters the optimized values 
become so high that the cost of the uranium necessary to make such a 
reactor critical would become prohibitively expensive. Also, a length 
to diameter ratio of approximately 5 and a power density of 9x10? W/m? 
would describe a highly inefficient reactor due to excessive neutron 


leakage through the core ends. Other important design factors (such as 
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shielding), which are not now contained in NUROC/SAC, indicate that the length 
should be only slightly larger than the diameter ele shoulda alcomb-egorea 

in the sensitivity data that wW is quite flexible (AW = +3x10’) and 

that smaller values may be chosen with little resulting loss in payload. 

This problem should serve as an example that one cannot randomly 
optimize variables or undertake a sensitivity analysis without first 
acquiring a knowledge of the system under consideration. With these 
facts in mind, a second optimization problem is solved in two-dimensions 


pmeem L/D fixed at 1.5. 


Maximize Payload: Free Parameters Q, W 


Starting Point; 


2000 MWt 


© 
Il 
~ 
Il 


6.0x10°? W/m? 


i 
I 
~ 
Ih 


“TABLE IV 


Sensitivity Data at the Optimum 









Ax for Specified Ay 


ndependent Simultaneous Variation 


“ip opt Vane Eom | A = 5151 A = 23630 














30,938Ke. Q=1823 MWt | AQ = +339. 340. a6 
=4.97x10° , | dv = 4.160 Ey £159 
W 
fia Calo alge) (x107) 


Analysis: 


The restriction on L/D (fixed at 1.5) results in a loss in optimal 
payload of only 2 percent. It is interesting to note that the optimal 
power and allowable deviation are almost identical to the three-variable 


problem. The power density however has a considerably lower optimal 
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value and AW is much smaller when L/D is fixed. 

The sensitivity data reveals that any power between 1500 and 2150 MWt 
may be used if w is held constant with only a small (1%) resulting loss 
in payload. In contrast the power density must remain close (within 47) 
to 5.0x10° W/M?. The fact that the eigenvectors are located close to the 
Q and wW axes implies that little more can be gained by varying Q and 


YW simultaneously. 


Maximize Payload: Free Parameters Q, L/D 

With the discovery that w and L/D cannot both be free and having 
fixed L/D at 1.5 it would now be advantageous to fix wW and optimize 
em QO and L/D. A value of 3.0x10° W/M is chosen for the power density 


.based on accepted values for existing nuclear rocket engines. 


Bearting Point: 
Q =x, = 2000 MWt 


L/D 2.0 


I 

Ma 
NO 

tl 


TABLE V 


Sensitivity Data at_ the Optimum 
Ax for opecitied ay 
Y AY = 1% X Independent Simultaneous Variation 
pe cue Variation X = G355men = 3e67- 







30,675 Kg. +397 AQ = 114 |AQ = 390 


ma is. AL/D = 0.79|AL/D = -.056 


a eg 


Analysis: 


When the power density is fixed both the power and the ratio of length 


to diameter have large acceptable variations. The power may vary from 1150 
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to 1950 with L/D and w fixed, and L/D may vary from 1.9 to 2.45 with 

Q and y fixed. This allows for considerable flexibility in designing 

an optimal engine as long as the power density remains relatively constant. 
The least sensitive direction of change is associated with the mini- 

mum eigenvalue. The associated eigenvector indicates that L/D may still 

vary from .9 to 2.45 and Q need not remain fixed, but may vary by +114 MWt 

as long as the direction of change from the optimal value coincides with 


changes in L/D. 


Mebit1cation of Results 

In order to verify that the predicted sensitivity of the optimized 
variables is accurate, the allowable deviations were substituted into 
NUROC/SAC and ESCAPE and the resulting payload was computed. A comparison 
was made to check if the payload remained within the predicted 1% of the 
maximum. The optimal values were first rounded-off to Q = 1550 +400 Mwt; 


> 


mme= 1.67 +.75. 
TABLE VI 


Comparison: 
Maximum payload minus 1% = 30,370 Kg. 


Independent Variation L/D Payload (Kg) 

Loe) De Oe ixed . 90 3027 
LAS B10) S70 

II. AQ, L/D fixed Q 
1150. MWt 30,295 
1950. MWt 30,420 

Simultaneous Variation Q L/D 
1435 90 S0ass 
66S Las B0R445 


aie 





The table of comparisons indicates that 1f 90) and )0/0) aresune:22-cq 
either independently or simultaneously the loss in payload is even less 
than the predicted 14. If decreased the deviations become slightly larger 
moan LA but are still highly accurate. This may not always be these ice. 

Lt must be realized that the Sensitivity analysis is accurdtcsen 
at the optimum. For functions which are very flat (i.e. insensitive to 
change) the predicted variations may be quite large, as indeed they are 
im the @xample given. When the variations are substituted, the function 
is no longer close to the optimum and results may vary considerably from 
Emose predicted. For this reason it 1s prudent to verify results especially 
at points far from the optimun. 

Another reason for verifying results is that there exist small in- 
herent errors in the optimization search, matrix inversion and eigen 
analysis subroutines. Care should also be taken in defining the conver- 
gence criterion because if.the computer is forced to make repeated searches 
near the optimum the values for the inverse Hessian matrix will be destroyed. 


Any resulting sensitivity analysis will be meaningless. 


Application of Penalty Factor 


In the previous two-dimensional optimization problem an optimal power 
of 1550 MWt and length to diameter ratio of 1.67 was established for a 
fixed power density of 3.0x10? M/W°. As was previously mentioned for 
equality constrained problems it is of interest to know the cost of the 
restriction on power density A. The following example employed the 


penalty factor method to optimize a three-variable problem with the constraint: 


ae 9 
G(x.) = X. BeOsae 0 


S/o 





peabtine Point. 
Q 
ED 
W 


Penalty Factor 


Results: 
Payload 


Q 


L/D 


The cost of keeping the power 


ae 8 
Per (10 W/m). 


=x, = 2000 MWt 


6.0 


lI 
*s 
I! 


=x, = 4,0x10° w/M? 


I 
Zo 
Il 


10,000 


= 30,670 Kg. 

= 1547 MWwt. 

= 17 

= 3.020x10° w/M? 
= G* - G = .020 


= 2Pe = 406. 


density constant is 406 Kg. of payload 


=e 
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IX. CONCLUSIONS AND RECOMMENDATIONS 


A sensitivity analysis at the optimum has been performed through 
the use of existing computer codes. By simply defining an objective 
function and any number of independent variables an optimized solution 
is found. In addition, it has been shown that the sensitivity of each 
of the variables at the optimum may be given in a useful format. With 
all others held constant, the range that each variable may assume 
before a specified degradation in performance occurs is given. Also, 
the length and direction of the axis of all ridges in the objective 
function are listed in order of decreasing magnitude. With this 
information at hand a complete knowledge of the sensitivity of all the 
input parameters is available to the user and decisions regarding 
optimal parameter choice can be made, taking into account the 
flexibility given by the pee knowledge. 

The methods employed were shown to be highly accurate when dealing 
with purely mathematical problems. The sensitivities of Rosenbrock's 
parabolic valley function were found with little difficulty. Engineering 
applications, on the other hand, require a great deal of insight into 
the problem before a sensitivity analysis may be started. Scaling the 
variables is important so that the sensitivity information is meaningful. 
The nuclear rocket example studied the engine performance related to 
power, power density and ratio of length to diameter of the nuclear core. 
A scaling problem existed because the power in watts and power density 
in watts/cubic meter were of the order of magnitude of iO while the 


length/diameter was approximately unity. Usually such a problem can be 





avoided by normalizing each variable by dividing it by the order of 
magnitude. This would be fine if the acceptable deviations occurred 
within the same number of significant figures. However, if one variable 
is much less sensitive to change than all the others, then the eigen- 
vector associated with the minimum eigenvalue will be dominated by that 
component and very little knowledge can be gained about the other 
variables in that direction. It would be convenient to modify the 
Computer program so that the sensitivities would be normalized instead 
of the variables. Of course, the user could always eliminate the 
problem by normalizing the variables the first time, and after observing 
the resulting deviations, normalize the sensitivities and optimize 
again. The optimized solution will be the same but the sensitivity 
data will be more accurate and meaningful. 

Another problem experienced when working with engineering problems 
was in defining a convergence criterion. For the Rosenbrock function a 
change in the variables of Ors produced significant changes in the 
objective function. When optimizing the nuclear rocker, however, the 
maximum payload was essentially determined when the normalized variables 
were accurate in the second decimal place. With the convergence criterion 
set at rome the program made over a were more iterations before 
stopping and found an optimum which was only about one kilogram more in 
payload. When the optimization search stays close to the optimum for 
many iterations, the inverse Hessian matrix is destroyed and any 
resulting sensitivity analysis has no meaning. A modification in the 
program is required so that when the optimum has essentially been 


determined, the inverse Hessian matrix can be stored and any subsequent 


refinement of the optimum will not affect the sensitivity data. 
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Equality constrained optimization has been discussed but not in 
detail. The penalty factor method used by Kelley was implemented 
to determine an approximation to what is usually referred to as the 
cost of the constraint. 

It would also be of interest to know how far in the plane of the 
constraint and normal to it, one may travel on the response surface 
before reaching a specified change in the optimum. This thesis has not 
considered such a sensitivity analysis. The techniques involved are 
similar, but more attention is needed in this area. 

Inequality constrained optimization problems are generally no more 
difficult to solve than equality constrained ones. For well behaved 
functions, the solution either does not violate the constraint or lies 
on the boundary and may be treated as an equality constraint. In 
either case, the techniques developed in this thesis can be applied. 
Care should be taken to try-a variety of starting points so that if the 
same solution is reached, the function can be assumed to be well 


behaved. 
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APPENDIX A 


This appendix describes the implementation of the variable metric op- 
timization method and matrix inversion and cigenvalue analysis computer 
programs for the purpose of sensitivity analysis. The input variables 
and their definitions are shown along with a flow diagram illustrating 
how the programs were modified to run successively and produce the desired 
sensitivity data. A typical output format is also given. 

The matrix inversion and eigenvalue analysis subroutine listings con- 
tain brief descriptions of the methods employed and their references. For 


a complete description of the computer coded variable metric method consult 


reference 6. 





BLOCK DIAGRAM OF CONTROLLING PROGRAM 


A-2 





> LARE 


2 ll 
LOAD 

MN aGLEP ao. 0 6. 

| NEXP, DELLAM 
CHANGE, JPOINT 
XJ), X()MIN, 
_X(J)MAX, DELX(J) ! 


— ] 
| 
[ee = 
{ 
| TEST FALSE | 
IN < 10 | Paes 
— i ) a oe 


X(J) > XMIN(J) 
X(J) < XMAX(J) | 
| TRUE 


$+ 





TPE NEIL 
MODEL (N,F,X) | 








INITIAL PRINTOUT 
NIMP, NGRD, NEVAL, 
i ey TST. 








SN 


OPTIMIZATION LOOP 


_FALSE XMIN(J)=X(J) 
XMAX (J) =X (J) 


CALL VARMET 


(NX, F, X, “MIN, XMAX, 
~ lap oy om 
‘eaok DELX Sine oc ee 


NEXP, DELLAM, CHANGE, 
JPOINT, JFLOW, H) 





a oe eee 


RETURN 


vA oe 
ve i \ 


TRUE 








yA a 
¢ JFLOW \ 
Ne 
Me 
FALSE 


‘PRINT OPTIMUM AND 
| CO-ORDINATES 
FB; .G)mie a 


a ee ee 


ae ee 


ee ee ee eee 


| | 
| TEST CONST RA TNS 


~——_ - 


TRUE 


y a . 


| PRINT INVERSE 
| HESSIAN MATRIX 
eel) 





_ CHANGE STORAGE 
| MODE 
HH(I) = H(I,J) 


AE 2 EAA DC 





—— 





v 


MATRIX 
INVERSION 





~ —— 


CALL MINV 


—_— 
¥ 


i (HH, N, D, L, 


oe 





ee ee 


| RETURN | | 


———_ 





_ PRINT HESSIAN 


| MATRIX 
| 


| HH(Z) 





oan 
j 


— ee 


CHANGE STORAGE | 


| 
| MODE 


——————————— eS 
| STORE DIAGONALS 
UNCOR(I) = H(I,I)_ 








L 


EIGENVALUE 
ANALYSIS 


< . 
G 





ye 
a 

Il 
, a 


4 


| CALL EIGEN 


?™ } 


M) j | CHE, R, N, MV) 


J 
me 


H(l,J) <= 


HH(I) 


eae 


— 


STORE EIGENVALUES 


| 


DECEIG(I) = 


H(I,J) = HH(I) | | 


—- 


} 


DY = 


ae. 2 eas 
a 








| ee: 











V. =k Vv 
ae 1 ne 
ones a 
‘ 
READY FOR OUTPUT 
PRINT 
OPTIMAL Y VECTOR OF OPTIMAL X 


ASSIGNED VARIATION 'DY' 
INDEPENDENT VARIATION =Dx 


SIMULTANEOUS VARIATION de> Vent 


BEALS 


RETURN 


CHANGE STORAGE MODE 


H(I,I) 


Fd 


ee — 


Rollo Live EOUA ELON S 


ee ae 
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Variablesietruc Program 


Input variables and definitions: 


N 
STEP 
ALPHA 
BETA 
NEXP 


DELLAM 


CHANGE 


JPOINT 


X(J) 
XMIN(J) 
XMAX (J) 


DECX(J) 


the number of independent variables 
initial step size for X components 


factor by which step size is increased (@ > 1) 


lv 


factor by which step size is decreased (0 < £6 < 1) 

limit on number of experiments along a vector line 
termination criterion for vector search expressed as a 
range fraction 

program termination criterion expressed as a range fraction 
control for type of numerical differentiation desired in 


gradient subroutine. 


i=) forward diiterence 
2 = backward difference 
pee Cental ditrercence 


vector of initial values for independent variables 

lower bound for each independent variable 

upper bound for each independent variable 

step size for each independent variable used in gradient 


subroutine for numerical differentiation. 


There are N + 2 data cards required as inputs for each program execution. 


The first is an identity card written in any format. The second card contains 


the first eight input variables listed above in the following format fields. 


(110, 3F10.4, 110, 2F10.4, 110) 


Each of the N remaining data cards contains the initial value, range 


and step size of the independent variables according to the following format 


statement. 


FORMAT (4F15.5) 





4 . ~ a 4 ~ iy * Awa «th Ta ~~ [=> ~~ - ro mo, LoS my Oy 
A typical output from the varigble metric optimizacion Sree se eee 


Zt should be noted that for the inverse Hessian matrix to be printedepro sem, 


the array designated H must be dimensioned exactly, i.e. if N = 2 then 


DIMENSION H (2,2). 





)BO 




























= 2 
‘BP Os C.10C00 
PHA = 1.00000 
q = C. 50000 
aAM = 4 01000 
ANGE = 0.00010 
OINT = 1 
a XMIN' ne XMAX DELY 
"-1.20000 -5.00000 5.00000 C.CO10C 
= 1.00000==——-5.00000-=— 5.¢0000 0.00100 
NIME NGRD NEVAL FUNCTION INDEPENDENT VARIABLES 
= 0.588002 01 =) 20C0 HMO. QmiC GC CimO 
==. 262008 00=S===- 0.511868 002-0. 3576 32-06 
Orb 0054s Te Or SIME 00 Oe 2a OO 
1 0.10833E-01 0.14259E 00 0. ao5ues =o 
we C. ante ene 059665220 70) m= Ald dose 
0 15 eae C8 6307 ae ree oe 
———— 0. 39 6G ee —0.20230E-06 Ti 218/22 ee 
TIMIZATION COMPLET. = SP TUNCTION EVIL DATION eee 
9 53 DOOR GO =0.202338-0iom e- 0. 32 cere 
HE 
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|  PRCGEAMME WAS EXECUTING LINE 65 IN ROUTINE M/PROG WHEN TERMINATION OCC 
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Matrix Inversion and Eigenvalue Analysis 

The example output of the variable metric computer code shown previously 
lists the optimum point and the inverse Hessian matrix at the optimum. 
Before the matrix inversion subroutine (MINV) may be implemented the upper 
Eriangular elements of the symmetric H matrix must be stored in 4 singly 
dimensioned array (HH). This was necessary because of the input format 
utilized in the calling sequence of MINV and was accomplished by inserting 
ae DO" loop in the main program. 

Once the matrix has been inverted, the elements of the Hessian are 
printed and the diagonal elements stored for calculation of sensitivity 
information. Implementation of the eigenvalue analysis subroutine 


(EIGEN) merely requires defining the constant MV, i.e. 


tt 


MV = 0 eigenvalues only 


MV 


1 eigenvalues and eigenvectors 

A dimension statement in the main program for sufficient storage space of 
the vectors L, M, (utilized by EIGEN) and R (utilized by MINV) is also 
necessary. All of the required information to construct the table of 
sensitivity data is now available in the main program. Insertion of the 
proper equations and output statements completes the modification for 
sensitivity analysis. 


A listing of MINV and EIGEN and an example output format follows. 
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